/////////////////////////////////////////////////////////////////////
//
// \author Jeongnim Kim  jnkim@ncsa.uiuc.edu
//
// Measuring abstraction penalty using TinyVector and Vector based on PETE
// src/PETE      (copy of PETE library http://www.acl.lanl.gov)
// src/OhmmsPETE (Ohmms implementation of a fixed vector and an array)
//////////////////////////////////////////////////////////////////////
#ifndef OHMMSCORE_VECTOR_ABSTRACTION_H
#define OHMMSCORE_VECTOR_ABSTRACTION_H

#include <config.h>
#include <OhmmsPETE/TinyVector.h>
#include <OhmmsPETE/OhmmsVector.h>
#include <Utilities/RandomGenerator.h>
#include <Utilities/Timer.h>
#include <Utilities/OhmmsInfo.h>

#if defined(ENABLE_FORTRAN) && defined(ADD_)
#define addmult addmult_
#endif
using namespace std;

extern void 
print_timing(const std::vector<double>& timingdata, int n, int niter, const std::string& helpmsg);

extern "C" 
{ // fortran call
  void addmult(double* a,const double* b,const double* c, const double *d,
	       int* ntot);
}

inline void 
caddmult(double* restrict a, const double* restrict b, 
	 const double* restrict c, const double* restrict d, int ntot) 
{
  for(int i=0; i<ntot; i++) {
     a[i] = b[i] + 1.3*c[i] + 1.7*d[i];
  } 
}


namespace  APPNAMESPACE
{
  template<class C, class RNG>
    class RandomVector{ };

  template<class T, unsigned D, class RNG>
    class 
    RandomVector<TinyVector<T,D>, RNG> {

      public:
        typedef TinyVector<T,D> Return_t;
        RandomVector(RNG& gen): RG(gen){ }
        inline Return_t operator()() {
          Return_t res(0);
          for(int i=0; i<D; i++) res[i] = 2*RG()-1;
          return res;
        }
        RNG& RG;
    };


  // specialized for TinyVector<T,D>
  template<class RNG>
    struct RandomVector<TinyVector<double,3>, RNG> {
      typedef TinyVector<double,3> Return_t;
      RandomVector(RNG& gen): RG(gen){ }

      inline Return_t operator()() {
        return Return_t(2*RG()-1, 2*RG()-1, 2*RG()-1);
      }
      RNG& RG;
    };


  template<class Pos_t, class PosArray_t, class ValArray_t>
    struct VectorAbstraction {

      RandomGenerator_t Random;
      std::string help;

      VectorAbstraction(const std::string& a):help(a){}
      ~VectorAbstraction(){}
      void run(int n, int niter) {

        if(niter == 0) niter=1;

        std::cout << "<!-- " << help << std::endl;
        std::cout << "Starting vector operations:" << n << " " << QMC_DIM <<"-dimensional elements" << std::endl;

        Random.init(0,1,-1); 

        RandomVector<Pos_t, RandomGenerator_t> rvgen(Random);

        double norm = 1.0/static_cast<double>(niter);
        APPNAMESPACE::Timer myclock;
        vector<double> timingdata(10);

        myclock.restart();
        PosArray_t A(n), B(n), C(n), D(n);
        timingdata[0] = myclock.elapsed();

        myclock.restart();
        for(int i=0; i<n; i++) {
          A[i]=0.0;
        }
        for(int i=0; i<n; i++) {
          B[i] = rvgen();
        }
        for(int i=0; i<n; i++) {
          C[i] = rvgen();
        }
        D = 1.1*B;

        timingdata[1] = myclock.elapsed();

        myclock.restart();
        for(int iter=0; iter<niter; iter++) {
          A = B + 0.13*C + 0.017*D;
        }
        timingdata[2] = myclock.elapsed()*norm;

        myclock.restart();
        for(int iter=0; iter<niter; iter++) {
          for(int i=0; i<n; i++)
            A[i] = B[i] + 0.13*C[i] + 0.017*D[i];
        }
        timingdata[3] = myclock.elapsed()*norm;

        if(n<=100) {
          for(int i=0; i<n; i++)   {
            std::cout << A[i] << " " << B[i] << " " << C[i] <<  " " << D[i] << std::endl;
          }
        }

        ValArray_t a(QMC_DIM*n), b(QMC_DIM*n), c(QMC_DIM*n), d(QMC_DIM*n);
        int ii=0;
        for(int i=0; i<n; i++) {
          for(int j=0; j<QMC_DIM; j++)
            b[ii++] = B[i][j];
        }

        ii=0;
        for(int i=0; i<n; i++) {
          for(int j=0; j<QMC_DIM; j++)
            c[ii++] = C[i][j];
        }

        ii=0;
        for(int i=0; i<n; i++) {
          for(int j=0; j<QMC_DIM; j++)
            d[ii++] = D[i][j];
        }

        myclock.restart();
        for(int iter=0; iter<niter; iter++) {
          a = b + 0.13*c + 0.017*d;
        }
        timingdata[4] = myclock.elapsed()*norm;

        int ntot = QMC_DIM*n;
        myclock.restart();
        for(int iter=0; iter<niter; iter++) {
          for(int i=0; i<ntot; i++) a[i] = b[i] + 0.13*c[i] + 0.017*d[i];
        }
        timingdata[5] = myclock.elapsed()*norm;

        myclock.restart();
        for(int iter=0; iter<niter; iter++) {
          caddmult(&(a[0]), &(b[0]), &(c[0]), &(d[0]), ntot);
        }
        timingdata[6] = myclock.elapsed()*norm;

#if defined(ENABLE_FORTRAN)
        myclock.restart();
        for(int iter=0; iter<niter; iter++) {
          addmult(&(a[0]), &(b[0]), &(c[0]), &(d[0]), &ntot);
        }
        timingdata[7] = myclock.elapsed()*norm;
#endif


        std::cout << "--> " << std::endl;
        print_timing(timingdata,n,niter,help);

      }

    };
}
#endif

/***************************************************************************
 * $RCSfile: vector_bench.h,v $   $Author: jnkim $
 * $Revision: 1.3 $   $Date: 2005/12/11 15:28:52 $
 * $Id: vector_bench.h,v 1.3 2005/12/11 15:28:52 jnkim Exp $
 ************************************************************************/
